Exploring the mechanism of Buyang Huanwu decoction in the treatment of lumbar disc herniation based on network pharmacology and molecular docking

Buyang Huanwu decoction (BYHWD), as one of the traditional Chinese medicine formulas, is widely used in the clinical treatment of lumbar disc herniation (LDH) with curative effect. It has the characteristics of multi-component, multi-target, and mutual synergy, but the mechanism of action is often unclear. It needs some research to explore the molecular mechanism of BYHWD in the treatment of LDH based on network pharmacology and molecular docking. Screen the active compounds of BYHWD and predict drug-related gene/protein targets, which could determine the specific target of BYHWD in the treatment of LDH. Construct the “Drugs-Compounds-Targets” network and search for the core targets. Use Gene Ontology functional enrichment analysis, Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis, and molecular docking verification to explore the possible molecular mechanism. Eighty-two effective compounds and 666 targets of BYHWD, 187 targets for LDH treatment, and 20 core candidate targets were excavated. A total of 3414 entries were identified by Gene Ontology enrichment analysis, 173 related signal pathways were identified by Kyoto Encyclopedia of Genes and Genomes enrichment analysis, and 5 core compounds were identified by molecular docking, which had a good affinity with core genes STAT3, JUN, AKT1, MAPK1, RELA, and PIK3CA. BYHWD may play the role of analgesic and improving function by synergistic anti-inflammatory and analgesic compounds, regulating cell metabolic differentiation, regulating immunity, and anticoagulation. BYHWD in the treatment of LDH may play a role in analgesia and improve function through multiple signaling pathways, including PI3K-Akt, mitogen-activated protein kinase, tumor necrosis factor, and interleukin-17. The PI3K-Akt signaling may be one of the key mechanisms.


Introduction
Lumbar disc herniation (LDH) is a clinical syndrome in which the annulus fibrosus of the lumbar disc is partially or completely ruptured due to various reasons, and nucleus pulposus tissue protrudes backward from the rupture opening, stimulating or compressing nerve root and cauda equina nerve. With the aging of population and the change of lifestyle, the incidence of LDH has increased significantly and tends to be younger. [1] LDH patients showed a variety of clinical manifestations such as low back pain, numbness, fatigue of lower limbs, etc. In severe cases, cauda equina nerve injury occurs, resulting in defecation disorder and even paralysis. The occurrence and progress of LDH have a serious impact on the quality of life, and it also brings serious economic pressure to patients' families due to the high disability rate. The treatment scheme of LDH includes conservative treatment and surgical treatment. Studies have shown that surgical treatment is more effective for patients with severe nerve root compression and inflammation symptoms, or without improvement after 3 to 6 months of non-surgical treatment. [2] However, conservative treatment is still the first choice for LDH. [3] Therefore, the development of complementary and alternative therapies for LDH is crucial.
Traditional Chinese medicine (TCM), one of the complementary therapies, has achieved significant effects in treating LDH and helped broaden the ideas of therapeutic approaches for LDH. Buyang Huanwu decoction (BYHWD), as one of the TCM formulas, including Huang-Qi, Dang-Gui, Tao-Ren, Hong-Hua, Chi-Shao, Chuan-Xiong, and Di-Long, is widely used in the clinical treatment of LDH with curative effect. [4,5] The TCM formula such as BYHWD has the characteristics of multi-component, multi-target, and mutual synergy, but the mechanism of action is often unclear. Network pharmacology screens synergistic multiple compounds from TCM formula in a high-throughput manner and explains the combination rules and network regulatory effects of TCM formulas. [6,7] We analyzed key targets and signal pathways of BYHWD in the treatment of LDH and explored its possible molecular mechanism based on network pharmacology, bioinformatics, and molecular docking (Fig. 1). and associated toolkits were used to analyze and complete visualization; Cytoscape 3.8.2 (National Institute of General Medical Sciences (NIGMS)) and related plug-ins are used to build visual networks; ChemBio 3D Ultra 14.0.0.117 (PerkinElmer) was used to convert the 2D structure to 3D form; PyMOL 2.4.0 (Schrodinger) was used to remove small-molecule ligands and water molecules and visualize them; Auto Dock Tools 1.5.6 (The Scripps Research Institute) software was used for hydrogenation and other processes; Vina 1.1.2 (Dr. Oleg Trott from The Scripps Research Institute) was used to process molecular docking and analyze binding activity; and GraphPad Prism 8.3.0 (GraphPad Software) was used to build heatmaps. This study was approved by the Ethics Committee of the First Affiliated Hospital of Zhejiang Chinese Medical University.

Main active ingredients and targets of BYHWD
Each Chinese herb included in BYHWD (Huangqi, Danggui, Taoren, Honghua, Chishao, Chuanxiong, and Dilong) was used as keywords for retrieval in the TCMSP database. The effective compounds of Chinese herbs were obtained from "Ingredients" and screened according to 2 criteria, "oral bioavailability ≥30%" and "drug-like index ≥0.18". [8,9] We collected potential protein targets for each drug in "Related Targets". For drugs that could not be retrieved in the database, their ingredients were obtained by searching the literature. We downloaded their structures from the PubChem database and entered them into the Swiss ADME platform. In the results of "pharmacokinetics", "GI absorption: high" was set as the precondition, and the standard that 2 or more other indicators display "yes" was used to screen compounds. [10,11] We subsequently obtained the compound corresponding acting target proteins by inputting the compounds in the Swiss Target Prediction platform and eliminated items with "probability" as "0". Finally, we filtered the items by "reviewed" in the UniProt database with the species selection "Homo sapiens (human)" to normalize the resulting protein targets.

Acquisition of disease targets for LDH
We selected "lumbar disc herniation" as the keyword and screened disease-related targets in GeneCards, OMIM, PharmGKB, DisGENET, and GAD databases. The Venn diagram was drawn by R 4.0.3 software and packages. We eliminated the duplicate items in Excel and obtained unique LDH-related targets.

Construction and analysis of visualized network
We used R software to match and map drug prediction targets to disease targets. We drew the Venn diagram of the "drugs-disease targets" intersection and extracted the common targetsthe key targets of BYHWD in treating LDH. Then we used Cytoscape to form a "Drugs-Compounds-Targets" visualization network.

Construction of protein-protein interaction network
We uploaded the key proteins obtained from network analysis to the online STRING 11.0 (ELIXIR infrastructure) database, setting the type as "homo sapiens", setting the confidence threshold as "highest confidence (0.900)", hiding the disconnected nodes, and kept other setting parameters unchanged to obtain protein-protein interaction (PPI) network.

Obtain key genes from PPI network
We imported the PPI network into Cytoscape and CytoNCA network topology analysis plug-in was used to calculate the betweenness centrality, closeness centrality, degree centrality, eigenvector centrality, local average connectivity, and network centrality. We screened out the nodes whose above indexes were greater than the median of all gene nodes of each index by R software. Then we created the sub-network, and key genes were calculated and screened again.

Enrichment analysis
We used R software "ClusterProfiler" package to perform Gene Ontology (GO) function enrichment analysis at the different biological process, cellular component, and molecular function levels of key targets. Then Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was carried out at key targets.

Verification of molecular docking
We downloaded the active ingredients' 2D structure of small molecule ligands from the PubChem database, and imported them into ChemBio3D Ultra, then calculated and set them as the 3D structure with minimum free energy. We downloaded 3D structure files of core target proteins from the PDB database and imported them into PyMOL software to remove small-molecule ligands and water molecules, then we used Auto Dock Tools software for hydrogenation. Finally, we performed molecular docking of ligand and receptor via Vina, analyzed binding activity, and constructed the heatmap with Prism.

Active ingredients and targets of BYHWD
According to the TCMSP database, the effective compounds of Astragalus Mongholicus, Angelica Sinensis, Radix Paeoniae Rubra, Peach Kernel, Safflower, and Ligusticum chuanxiong were 119, 66, 189, 87, 189, and 125, respectively. Then, based on oral bioavailability ≥30% and drug-like ≥0.18, 103 compounds were screened out, including 20 of Astragalus Mongholicus, 2 of Angelica Sinensis, 29 of Radix Paeoniae Rubra, 23 of Peach Kernel, 22 of Safflower, and 7 of Ligusticum chuanxiong. Eighty-six compounds of "Earthworm" were obtained through literature reports, [12] and 20 effective compounds were screened on the SwissADME (Swiss Institute of Bioinformatics) platform. Finally, we obtained a total of 82 effective ingredients by combining all the active ingredients and removing repetitive items. A total of 3107 targets were corresponding to the effective ingredients of these 7 drugs. The protein targets were standardized, and 666 unique protein targets were obtained after removing unrecognized and duplicate items.

Disease targets of LDH
A number of 396, 17, 1658, 1, and 3 disease targets of LDH were obtained in the databases of Gene Cards, OMIM, Pharm GKB, DisGeNET, and GAD, respectively, and 1953 LDHrelated targets were obtained after merging and de-duplication ( Fig. 2A). After mapping the drug targets to the LDH targets, 187 common targets were obtained, such as PTGS2, NOS2, IL6, and MAPK1 (Fig. 2B).

Visualization network
We used Cytoscape to draw the relationship between the above active ingredients and targets. Two hundred fifty-one nodes and 959 relationship edges were obtained, and a visual network diagram of "Drugs-Compounds-Targets" was constructed. Compound nodes were presented as pie charts and different colors, and different colors represent that the compounds come from multiple Chinese herbs and drugs (Fig. 3).

PPI network
We uploaded the screened common protein targets to the STRING 11.0 platform and set the type as "homo sapiens", then the confidence threshold was set as "the highest confidence (0.900)" according to the complexity of the network. We hid the disconnected nodes and kept other setting parameters unchanged, and ultimately constructed the PPI network diagram. The PPI network contains 187 proteins, 995 action edges, and 324 expected edges with enrichment P value <1.0e −16 (Fig. 4A).

Find core genes of PPI network
The above data files were imported into Cytoscape to generate a protein interaction network diagram, and the betweenness centrality, closeness centrality, degree centrality, eigenvector centrality, local average connectivity, and network centrality were calculated by CytoNCA network topology analysis plug-in, and finally, 175 nodes and 995 interaction edges were obtained. The nodes with the above indexes greater than the median of all gene nodes of each index were screened out, and a total of 52 nodes and 402 edges were obtained. Then the matrix was imported into Cytoscape to extract sub-networks, and 20 key nodes and 224 edges were calculated and screened again. Finally, the core network was extracted, and the nodes were key genes. The corresponding indexes were calculated in descending order of degree (Table 1). Therefore, it was predicted that STAT3, JUN, Medicine AKT1, MAPK1, and RELA in BYHWD were the key genes for treating LDH (Fig. 4B).

GO function enrichment analysis
GO enrichment analysis identified 3414 items of BYHWD in treating LDH and listed the top enrichment analysis results respectively (q value < 0.05). Its biological processes mainly include epithelial cell proliferation, response to lipopolysaccharide, response to bacterial molecules, gland development, and so on. Cell compounds mainly include membrane raft, membrane microdomain, membrane region, transcription regulatory complex, focal adhesion, and so on. Molecular functions mainly include DNA and RNA transcription factor binding, protease activity, cytokine receptor binding, nuclear receptor activity, transcription factor activity, etc (Fig. 5A).

KEGG signaling pathway enrichment analysis
The results of KEGG enrichment analysis showed that the key targets of BYHWD acting on LDH were mainly related to 173 signaling pathways such as PI3K-Akt signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway, hepatitis B, tumor necrosis factor (TNF) signaling pathway, interleukin (IL)-17 signaling pathway, etc. Thirty pathways with a large number of enriched genes were screened out, of which the PI3K-Akt signaling pathway was the most significantly enriched (Fig. 5B and C).

Molecular docking results
It is generally believed that the lower the energy is spent, the more likely the interaction will occur when the conformation of the ligand bound to the receptor is stabilized. The binding energy is less than -4.25 kcal mol −1 , indicating that the ligand has a certain binding activity with the receptor. The binding activity is better when it is less than -5.0 kcal mol −1 , and strong when it is less than -7.0 kcal mol −1 . [13] The top 8 active ingredients of the "Drugs-Compounds-Targets" network, including quercetin, kaempferol, beta-sitosterol, (+/-)-18-hydroxy-5Z, 8Z, 11Z, 14Z, 16E-eicosapentaenoic acid, luteolin, tridecanoic acid, arachidonic acid, and baicalein, were individually docked with the top 6 core protein targets, including STAT3, JUN, AKT1, MAPK1, RELA, and PIK3CA, and the binding energy threshold was set as less than -5.0 kcal mol −1 . The results showed that quercetin, kaempferol, beta-sitosterol, luteolin, and baicalein were successfully docked to the corresponding proteins with good stability. Among them, the strongest binding was between beta-sitosterol and AKT1 (-10.8 kcal mol −1 ), and except the binding energy between beta-sitosterol and STAT3 was -6.5 kcal mol −1 , all others had strong binding activity. This indicated that these compounds played a role through multiple targets. The binding activity of tridecanoic acid was the weakest in the overall level. Its binding activity with AKT1 was -5.6 kcal mol −1 , which was stable, and the binding activities with other genes were greater than -5.0 kcal mol −1 , which were relatively poor. All the compounds were successfully docked with AKT1 and had good binding activity (Fig. 6A). We selected the top 3 compounds including quercetin, kaempferol, and beta-sitosterol to visualize the docking results with AKT1 and output a 3-dimensional diagram of molecular docking ( Fig. 6B-E).

Discussion
At present, the pathogenesis of LDH is gradually being recognized which includes mechanical compression of the nerve root, inflammatory stimulation of nerve root, apoptosis, and autoimmune mechanism of intervertebral disc. [14][15][16][17][18] However, with the aggravation of the aging of the social population, the occurrence of LDH has a trend of younger age. Therefore, early detection and intervention treatment should be actively carried out. At present, conservative treatment mainly includes oral drug therapy, percutaneous block therapy, physical therapy, and exercise therapy. [19,20] Among them, oral drugs such as non-steroidal anti-inflammatory drugs have better curative effects, but also have significant toxic and side effects such as gastrointestinal reactions and nephrotoxicity. The long-term application can lead to gastric mucosal erosion, gastric hemorrhage, peptic ulcer formation, and even perforation death. Recently, the advantages of Chinese herb are constantly being reflected, but it has the characteristics of multi-compounds, multi-targets, and multi-pathways, and all compounds work synergistically, so it is difficult to clarify the material basis and specific mechanism of action.
We constructed a "Drugs-Compounds-Targets" network and the above results showed that the effective components can improve LDH symptoms and the progress of intervention, including quercetin, kaempferol, beta-sitosterol, luteolin, and baicalein. Many recent studies have found that immune cells and immune factors produced by them were an important immunological mechanism of nerve injury in LDH. [21][22][23] The expressions of inflammatory cytokine IL-1β, IL-6, and TNF-α were markedly increased in animal models of LDH, whereas these inflammatory cytokines have significant hyperalgesic effects. Quercetin reduces the production of inflammatory cytokines by inhibiting the MAPK/toll like receptor 4 signaling pathway and reduces the polarization of helper T cells 17 to exert anti-inflammatory effects. [24] Nowadays, studies have found that kaempferol inhibits oxidative stress and neuroinflammation. [25] The anti-inflammatory mechanisms of beta-sitosterol may differ from those of NSAIDs. It promotes cellular calcium uptake, inhibits the activity of myeloperoxidase and adenosine deaminase, and inhibits IL-1β, TNF-α levels. [26,27] Luteolin is widely distributed in nature and numerous Chinese medicinal herbs contain this component, such as Lonicera japonica, Prunella vulgaris, and guar vegetables, which have anti-inflammatory, expectorant antitussive, antitumor, and antiviral pharmacological activities. [28] Baicalein significantly suppressed the release of NO, IL-6, IL-1β, and TNF-α, thereby inhibiting inflammation occurrence. [29] It has been demonstrated that baicalein can be converted to baicalin in vivo, which inhibits the production of inflammatory factor TNF-α, IL-1β, IL-6, IL-17, and matrix metalloproteinase-9, and regulates NF-κB signaling pathway to exert anti-inflammatory effects. [30] In this study, 82 active compounds of BYHWD were excavated. A total of 3107 targets were corresponding to the   effective components of the 7 drugs, and 666 unique protein targets were obtained after standardization by the Uniprot database. There were 187 possible targets for the treatment of LDH, of which 20 were core targets such as STAT3, JUN, AKT1, MAPK1, RELA, etc. The results of GO enrichment analysis showed that BYHWD participated in signal transduction, cell metabolism, DNA and RNA transcription, and other processes. As a signal transducer and activator of transcription, the activation of STAT3 is considered to be closely related to the development of chronic pain and participates in the pathway regulating the occurrence of intervertebral disc degeneration. [31,32] The JUN gene is mainly involved in signal pathways such as apoptosis inhibition and neuronal degeneration. [33,34] AKT1 is a core factor in the PI3K-Akt signaling pathway and plays an important role in cell growth, metabolic regulation, and cancer. [35] MAPK1 is mainly involved in a variety of cellular processes, such as proliferation, differentiation, transcriptional regulation, and development. [36] RELA is an important member of the NF-κB family. Its post-translational modification can accurately regulate the transcription activity of NF-κB and plays an important role in regulating important life activities such as inflammation, tumor, metabolism, and immune response. [37] Therefore, it was speculated that BYHWD mainly alleviates the symptoms of patients by participating in cell metabolism, inflammatory reaction, transcriptional regulation, and immune response. One hundred seventy-three signaling pathways were found by KEGG enrichment analysis. The pathways with higher gene ratio include PI3K-Akt signaling pathway, MAPK signaling pathway, virus infection-related pathway, tumor-related pathway, diabetic complication-related pathway, TNF signaling pathway, IL-17 signaling pathway, etc. Among them, the PI3K-Akt signaling pathway was the most significant enrichment. LDH has no significant correlation with tumor-related pathways and virus infection-related pathways, so we would not pay attention to them for the time being. IL-17A promotes intervertebral disc degeneration by inhibiting autophagy through activation of the PI3K/Akt/Bcl-2 signaling pathway. [38] The MAPK pathway participates in central sensitization and nociceptive specific signaling in dorsal horn neurons. [39] This study focuses on the PI3K-Akt signaling pathway. The PI3K-Akt signaling pathway is an important node in mammalian cells that controls cell growth, migration, proliferation, and metabolism, and plays an important role in regulating T cell development, function, and stability. [40] The gene expression results studied by Tian et al show that hypoxia and nutritional deficiency can increase the relative expression of PI3K and Akt genes and inhibit the expression of functional genes. [41] PI3K-Akt signaling pathway has a protective effect on mesenchymal stem cells derived from human nucleus pulposus against hypoxia and nutritional deficiency. When the PI3K-Akt pathway is blocked, it will aggravate intervertebral disc degeneration caused by inflammation, and inhibit the PI3K-Akt signal to promote osteogenic differentiation of mesenchymal stem cells in osteoporosis rat model. [42,43] Moreover, the PI3K-Akt-mTOR signaling pathway is important for the normal metabolism of joint tissue. [44] Inhibition of PI3K-Akt-mTOR signaling pathway leads to  inhibition of fibroblast proliferation reduction, stimulation of apoptosis and autophagy, and provides new ideas for reducing epidural fibrosis caused by surgery. [45] Xu et al [46] found that LDH significantly increases thermo-mechanical abnormal pain and leads to degeneration of dorsal root ganglion (DRG) by inactivating the L-PGS/PI3K/Akt pathway. Upregulation of Nav1.6 can increase neuronal excitability and participate in neuropathic pain of DRG. By activating the TNF-/STAT3 pathway, histone H4 is superacetylated in the STAT3-mediated DRG Scn8a promoter region, resulting in L5-VRT-induced neuropathic pain. [47] However, there are also deficiencies. First, the information of drug ingredients in databases and literature may still not be exhaustive, and insect drugs have a different research system than herbs. Then again, this study was only corroborated by relevant literature. Therefore, in a follow-up study, our group intends to intervene biological samples with full prescription drug delivery and decoupled drug delivery, and carry out high-throughput assays such as transcriptomics and proteomics to obtain the effect target profiles of BYHWD, to provide a high-quality data resource for more comprehensive annotation of the action mechanism of BYHWD.

Conclusions
In this study, we explored the mechanism, core targets, and signaling pathways of BYHWD in the treatment of LDH through network pharmacology and bioinformatics analysis. It was speculated that BYHWD plays a role in the treatment of LDH through multiple signaling pathways such as PI3K-Akt, MAPK, IL-17, and TNF, and PI3K-Akt signaling pathway may be one of the key mechanisms. These were conducive to optimize the experimental design, provide a reference for new drug development, and save costs to a certain extent. We also hope to promote BYHWD as a supplement and alternative treatment for LDH.

Author contributions
GY wrote the paper. ZHJ, WXJ, ZSX, and LSJ performed the literature search and revised the paper. TPJ and LSJ guided the writing of the paper and reviewed the manuscript. All authors critically reviewed progressive drafts of the manuscript and approved the final version.